use dataset_county, clear

egen cid = group(statefips countyfips)

keep D199* D20* cid
reshape i cid
reshape j year
reshape xi
reshape xij D
reshape long

collapse (sum) D, by(year)

keep if year < 2017

replace D = D / 1000

#delimit;

gr tw
	(line D year, col(black) lwid(thick))
	,
		ylab(, angle(horiz))
		plotregion(style(none))
		xtitle("Year")
		ytitle("Number of wells (1,000s)")
		;
		
#delimit cr

gr export "_output/figure1.pdf", replace

